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ABSTRACT 

We present the latest development of the disk gravitational instability and fragmentation model, 
originally introduced by us to explain episodic accretion bursts in the early stages of star forma¬ 
tion. Using our numerical hydrodynamics model with improved disk thermal balance and star-disk 
interaction, we computed the evolution of protostellar disks formed from the gravitational collapse 
of prestellar cores. In agreement with our previous studies, we find that cores of higher initial mass 
and angular momentum produce disks that are more favourable to gravitational instability and frag¬ 
mentation, while a higher background irradiation and magnetic fields moderate the disk tendency to 
fragment. The protostellar accretion in our models is time-variable, thanks to the nonlinear inter¬ 
action between different spiral modes in the gravitationally unstable disk, and can undergo episodic 
bursts when fragments migrate onto the star owing to the gravitational interaction with other frag¬ 
ments or spiral arms. Most bursts occur in the partly embedded Class I phase, with a smaller fraction 
taking place in the deeply embedded Class 0 phase and a few possible bursts in the optically visible 
Class II phase. The average burst duration and mean luminosity are found to be in good agreement 
with those inferred from observations of FU-Orionis-type eruptions. The model predicts the existence 
of two types of bursts: the isolated ones, showing well-defined luminosity peaks separated with pro¬ 
longed periods (~ 10'’’ yr) of quiescent accretion, and clustered ones, demonstrating several bursts 
occurring one after another during just a few hundred years. Finally, we estimate that 40%-70% of 
the star-forming cores can display bursts after forming a star-disk system. 

Subject headings: accretion, accretion disks - hydrodynamics - instabilities ~ ISM: clouds - stars: 
formation 


1. INTRODUCTION 

Low-mass stars form as a result of the gravitational 
collapse of dense gaseous cores. Standard models of 
core collapse predict that the accretion rate onto a form¬ 
ing pr otostar is proportional to the cube of t he sound 
speed (|Larsonl [19691 [PenstonI [To^ IShd [19771 1. When 
a finite size of the core is taken into account, numer¬ 
ical simulations of gravitationally unstable, spherically 
symmetric cores indicate that accretion is tapering off 
with time in the late evo lution (jFoster fc Chevali^l 199,11 
IVorobvov fc Basull2005ll . However, when even a modest 
degree of rotation is present initially in the cor e as sug¬ 
gested by observations (e.g. ICaselli et al.ll200^ . simple 
arguments based on the centrifugal radius and sophisti¬ 
cated numerical hydrodynamics simulations both demon¬ 
strate that most of the core mass does not fall directly 
onto the protostar but rather lands onto an accretion 
disk formed from conservation of angular momentum of 
the core. 

It has recently become evident that the mass infall 
rate onto the disk Minfaii at radial scales on the order 
of 1000 AU and the mass accretion rate onto the star 
M may be significantly different thanks to the compli¬ 
cated interplay of various physical mechanisms of mass 
and angular momentum transport operating in the disk. 
For instance, Minfaii in isolated core models gradually 
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declines with time from a few x 10 ® Mq yr ^ to a neg¬ 
ligible value by the end of the embedded phase, but M 
in gravitationally unstable disks can be highly variable 
(|Vorobvovll2009l : lRice et al.ll201Cll) . sometimes exhibiting 
episodic accretion bursts > 10“^ Mq yr“^ caused by 
disk gravitational fragmentatio n and migration of the 
fragments on to the protostar (|Vorobvov fc Basiil 120051 
l2006l l2010af) . The formation of giant planets in the 
disk can also significantly alter the character of accre¬ 
tion, creating various patterns of variabi lity and bursts 
exceeding in magnitude 10~^ M g^ yr“^ (|Machida et al.l 
120111 : iNavakshin fc Lodatcill2012f) . A combination of the 
magneto-rotational and thermal instabilities in the in¬ 
ner several AU and gravitational instability further out 
in the disk (or layered accretion) was shown to produce 
accretion bursts typical for FU-Orion is-type eruptions 
(|Armitage et al.l l2001l : IZhu et akl l2009tl . To complicate 
the things further, Minfaii may experience significant vari¬ 
ations when the chaotic and turbulent nature of clustered 
star formation is taken int o consideration le.g. lBate et al.1 
IMfilPadoan et al.l[20ll . 

Observations support the growing evidence that accre¬ 
tion onto low-mass protostars is at least partly variable. 
An ever growing number of FU-Orion i s-type and EX- 
Lupi-like eruptive stars (|Audard et al.l l2014ll does not 
ht into the standard models of spherical core collapse. 
The mean/median luminosity of protostars in young star¬ 
forming regions appears to be lower by about an order of 
mag nitude than that predicted from the stan dard models 
(e.g. iKenvon et al.l Il99fll : lEvans et al.l l2009t) . Accretion 
rates gradually declining in time and showing episodic 
bursts were shown to resolve this ’’luminosity problem” 
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(|Dunham fc VorobvovI [2M^ . Monitoring of the accre¬ 
tion variability suggest that about half of all protostars 
show up to 50% variatio ns in M over timescales less than 
2 yr (|Billot et al.ll2012h . Surveys of young stars in star¬ 
forming regions and in the Galactic plane revealed that 
about 0 . 1 % of objects show a l uminosity increase by more 
than a factor of 2.5 over 5 yr (IScholz e t al.ll201 3H. includ- 
ing recent flares in VSX J205126.1 ( Kosnal et al.l[20TiIl 
and V1647 Ori (I Abraham et al.ir2004D . The knotted mor¬ 
phology of jets seen in some protostellar systems suggests 
an underlying variability in the mass accretion, although 
the combination of jet velocities and spacing between the 
knots often suggest shorter periods of episodidty th an 
predicted for FU-Orionis-type stars (jArce et al.ll2013ll . 

To summarize, accretion onto young stars seems to 
exhibit a variety of patterns with time variations of 
different amplitude a nd d uration and, as noted in 
iDunham fc VorobvovI (|2012ll . is better termed as vari¬ 
able accretion with episodic bursts. This newly emerg¬ 
ing paradigm is beginning to supersede the classical Shu- 
Larson-Penston steady accretion models and may have 
important consequences for the evolution of stars and 
planets. For instance, variable accretion can help to 
explain the luminosity spread of y oung clusters with¬ 
out in voking a significant age spread (|Baraffe et al.ll200M 
l2012t) . In addition, quiescent periods between accretion 
bursts can promote disk fragm entation and giant planet 
formation (|Stamatellosl 120111) . Finally, variable accre¬ 
tion with episodic bursts is expected to have a signifi¬ 
cant impact on the disk and envelope chemis try and on 
the composition of ices in protostellar disks (|Leell200^ 
IVisser fc Berginl[2012t iKim et ani2012tl 

In this paper, we revisit the disk instability and frag¬ 
mentation model for episodic accretion and luminosity 
bursts, originally d evelo ped by us in a series of papers 
(IVorobvov fo Basnl[2nnaf2(Ml2ninall . using an improved 
numerical hydrodynamics code which takes into account 
a better disk thermal physics, improved dust opacities, 
and an accurate calculation of the stellar photospheric 
and accretion luminosities using a stellar evolution code 
that takes stellar accretion into account. The latter 
update allows us to calculate the burst statistics and 
perform direct comparison with observations, and also 
make prediction regarding the expected fraction of star¬ 
forming cores that can display bursts after forming a 
star-disk system. The paper is organized as follows. A 
brief description of the numerical model and recent up¬ 
dates are presented in Section [21 The main results are 
described in Section [3] The characteristics of the bursts 
obtained in the framework of our model are reviewed 
in Section 01 The time evolution of individual bursts is 
considered in Section I The expected fraction of star¬ 
forming cores than can exhibit bursts after forming a 
star-disk systems is calculated in SectionjTland main con¬ 
clusions are summarized in Section 01 

2. MODEL DESCRIPTION 


Our numerical mode l is described in detail in 
IVorobvov fc BasiJ (|2010ah and is briefly reviewed below 
with the emphasis on several recent updates. We start 
our numerical simulations from the gravitational collapse 
of a starless cloud core, continue into the embedded phase 
of star formation, during which a star, disk, and envelope 
are formed, and terminate our simulations when the age 


of the star becomes older than 1.0 Myr. Such long inte¬ 
gration times are made possible by the use of the thin- 
di sk approximation, the justi fication of which is provided 
in IVorobvov fc Basul (l2010all . The protostellar disk oc¬ 
cupies the inner part of the numerical polar grid and is 
exposed to intense mass loading from the infalling enve¬ 
lope. 

To avoid too small time steps, we introduce a “sink 
cell” at Tsc = 6.0 AU and impose a free inflow inner 
boundary condition and a free outflow outer boundary 
condition so that the matter is allowed to flow out of the 
computational domain but is prevented from flowing in. 
The sink cell is dynamically inactive; it contributes only 
to the total gravitational potential and secures a smooth 
behaviour of the gravity force down to the stellar surface. 
During the early stages of the core collapse, we monitor 
the gas surface density in the sink cell and when its value 
exceeds a critical value for the transition from isothermal 
to adiabatic evolution, we introduce a central point-mass 
object. In the subsequent evolution, 90% of the gas that 
crosses the inner boundary is assumed to land on the 
central object. The other 10% of the accreted gas is 
assumed to be carried away with protostellar jets. 

2 . 1 . Main equations 

The basic equations of mass, momentum, and energy 
transport in the thin-disk limit are 

f)y 

— =-V,.(E«p), (1) 


dt 
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( 3 ) 

where subscripts p and p' refers to the planar compo¬ 
nents (r, (j)) in polar coordinates, E is the mass surface 
density, e is the internal energy per surface area, V is 
the vertically integrated gas pressure calculated via the 
ideal equation of state as 7^ = (7 — l)e, Z is the radially 
and azimuthally varying vertical scale height determined 
in each computational cell using an assumption of local 
hydrostatic equilibrium, Vp = VrV + v^cj) is the veloc¬ 
ity in the disk plane, and Vp = rdjdr + (f)r~^d/d(j) is 
the gradient along the planar coordinates of the disk. 
The gravitational acceleration in the disk plane, Qp = 

grV + gci)<p, takes into account self-gravity of the disk, 
found by solving for the Poisson integral (see details in 
IVorobvov fc Basull2010all . and the gravity of the central 
protostar when formed. Turbulent viscosity is taken into 
account via the viscous stress tensor 11 , t h e expr ession 
for which is provided in IVorobvov fc Basul (|2010ar) . We 
parameterize the magnitude of kinematic viscosity v us¬ 
ing the a-prescription with a spatially and temporally 
uniform a. 

Two l.h.s. terms in Equation 0 represent mag¬ 
netic pressure and tension in the thin-disk approxima¬ 
tion, where Bz is the vertically uniform magnetic field 
in the disk and Bp = B^r -|- B'^ij) are the planar 
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components of the magnetic held at the top surface of 
the disk. In the hux-freezing approximation adopted 
in this work the vertical magnetic held component in 
the disk can be determined from the r elation = 
2'iTG^^‘^T,j (|Nakano fc Nakamuralll97^ . where is 
the spatially uniform mass-to-hux ratio. The planar 
components of the magnetic held are directly related to 
the planar components of gravitational acceleration 
through the following relation B^Bp /{2'k) = —g^jg^ 
Isee I Vorobyov h Basull 200 ^ for more details). 

The radiative cooling A in equation Q is determined 
using the diffusion approximation of the vertical radia¬ 
tion trans port in a one-zone rnodel of the vertical disk 
structure (|johnson h Gamini3l2003ll 


A _ TT ^ T'4 ‘ 


(4) 


where r is the optical depth to the disk midplane, cr is 
the Stefan-Boltzmann constant, T^p = Vg/KS is the 
midplane temperature of gas^, g = 2.33 is the mean 
molecular weight, R is the universal gas constant, and 
y, = 2 -|- 20 tan“^(r)/( 37 r) is a function that secures a 
correct transition between the optically thick and opti¬ 
cally thin regimes. The heating function is expressed as 

r = (5) 

L-\-T 


where Tjir is the irradiation temperature at the disk sur¬ 
face determined by the stellar and background black- 
body irradiation as 

( 6 ) 

where Tbg is the uniform background temperature (in our 
model set to the initial temperature of the natal cloud 
core) and is the radiation flux (energy per unit 

time per unit surface area) absorbed by the disk surface 
at radial distance r from the central star. The latter 
quantity is calculated as 

^irr(?') = COS7irr, (7) 


where 7 irr is the incidence angle of radiation arriving at 
the disk surface (with respect to the normal) at radial 
distance r. The stellar luminosity L* is the sum of the 
accretion luminosity L*,accr = (1 — e)GM*M/2i?* aris¬ 
ing from the gravitational energy of accreted gas and the 
photospheric luminosity T*,ph due to gravitational com¬ 
pression and deuterium burning in the stellar interior. 
The stellar mass M^, and accretion rate onto the star M 
are determined using the amount of gas passing through 
the sink cell, while the stellar radius i?* is returned by 
a stellar evolution code (see Section 12.21 for details and 
the definition of e). Equations (P)-® are solved in the 
polar coordinates on a numerical grid with 512 x 512 grid 
zones. The solution proc edure is described in detail in 
IVorobvov h Basul (l2010all . 


2 . 2 . Recent updates 


^ This definitio n of the midplan e temperature is accurate within 
a factor of unity llZhu et al.ll2012h 


In this study, several important updates have been im¬ 
plemented to the numerical code as c ompared to the ear¬ 
lier work of IVorobvov Sz BasrJ (l2010all . Fi rst, we have im¬ 
plem ented newer Semenov dust opacities (ISemenov et al. 
2QQ^ instead of older Bell & Lin opacities ( Bell fc Lin 


1994fl . The Semenov opacities are somewhat higher than 


those of Bell & Lin in the temperature range typical for 
protostellar disks, which results in a somewhat higher gas 
temperature. Second, we considered a stiffer equation 
of state taking into account the fact that the rotational 
and vibrational degrees of freedom of molecular hydrogen 
are ex cited only above 100 K le.g. lMasnnaga fc Inntsukal 
1200011 . As a result, the ratio of specific heats takes the 
following from 


( 5/3, if Tg < 100 K , 

7 = < 7/5, if 100 K < Tg < 2000 K, ( 8 ) 
\ 1.1, ifTg> 2000 K. 

In IVorobvov fc Basul (j2010all . 7 was set to 7/5 for Tg < 
2000 K. The net result is an overall moderate increase in 
the disk temperature at distances > 10 AU. 

The aforementioned updates enable a better calcula¬ 
tion of the thermal balance in the disk and hence a more 
accurate study of the gravitational instability and frag¬ 
mentation. We note that the net effect of these updates 
is an increase in the disk temperature, making disk frag¬ 
mentation more difficult. 

Finally, we have improved on the method in which 
the parameters of the ce ntral star, are calculated. In 
IVorobvov fc BasiJ (|2010al l. pre-main-sequence tracks for 
the non-accreting low-mass s tars and brown dwarfs of 
iD’Antona &: Mazzitellil (Il994ll were employed to calcu¬ 
late the stellar photospheric luminosity and radius. In 
this work, the properties of the forming protostar are 
calculated using a stellar evoluti on code descr i bed i n 
iBaraffe k, Chabrieil (|2010t) . As in iBaraffe et al.l (|2012f) , 
we assume that the fraction e of the accretion energy 
GM*M/(2i?*) is absorbed by the protostar, while the 
fraction (1 — e) is radiated away and contributes to the 
accretion luminosity of the star T*_accr- Despite many 
efforts, the exact value of e in low-mass star formation is 
not known. In the present c alculations, we ado pt a so- 
called ’’hybrid” scheme ((see IBaraffe et al.ll 2012 l for de¬ 
tail) with e = 0 when accretion rates remain smaller than 
a critical value Mcr = 10 “® Mq yr“^, and e = 0.2 when 
M > Mcr- 

The stellar evolution code is coupled with the main hy- 
drodynamical code in real time. The input parameter to 
the stellar evolution code provided by disk modeling is 
the mass accretion rate onto the star M. The output of 
the stellar evolution code are the stellar radius i?* and 
the photospheric luminosity T,_ph, which are employed 
by the disk hydrodynamics simulations to calculate the 
total stellar luminosity and the radiation flux reaching 
the disk surface. Due to heavy computational load the 
stellar evolution code is invoked to update the properties 
of the protostar only every 20 yr, while the hydrodynam- 
ical time step may be as small as a few weeks and the 
entire duration of numerical simulations may exceed 1.0 
Myr. 

The coupling of disk modeling with the stellar evo- 
lution code i s essential . As was demonstrated in 
IBaraffe et al.l (l2009l[2M^ . the stellar properties derived 
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TABLE 1 

Model parameters 


Model 

^core 

(Mq) 

/3 

(%) 

Tinit 

(K) 

£^0 

(km pc~^) 

ro 

(AU) 

So 

(g cm-2) 

-Rout 

(pc) 

a 


1 

1.1 

0.88 

10 

1.4 

2400 

5.2 X 10“^ 

0.07 

5 X 10-1* 

0 

2 

1.5 

0.88 

10 

1.0 

3400 

3.7 X 10-2 

0.1 

5 X 10-3 

0 

3 

0.31 

0.88 

10 

3.3 

685 

1.8 X 10-1 

0.03 

5 X 10-3 

0 

4 

1.1 

0.14 

10 

0.57 

2400 

5.2 X 10-2 

0.07 

5 X 10-3 

0 

5 

1.1 

0.88 

25 

5.7 

960 

3.2 X 10-1 

0.028 

5 X 10-3 

0 

6 

1.1 

0.88 

10 

1.4 

2400 

5.2 X 10-2 

0.07 

3 X 10-2 

0 

7 

1.1 

0.88 

10 

1.4 

2400 

5.2 X 10-2 

0.07 

5 X 10-3 

3.33 


from evolution mode ls that do not take stellar accr etion 
into account (such as iD’Antona fc Mazzitellil (jl994l) data 
used in our previous work) can significantly differ from 
those derived from accreting models. For instance, the 
stellar photospheric luminosities may differ by a factor 
of up to 10 (see fig. 5 in Baraffe et al. 2012) and stellar 
radii® by a factor of several (see fig. 2 in Baraffe et al. 
2009), depending on the stellar mass and the fraction of 
accretion energy absorbed by the star e. Therefore, em¬ 
ploying accreting stellar evolution models is crucial for an 
accurate comparison of burst characteristics derived from 
numerical modeling with those measured in young star¬ 
forming regions (Sections |3][7]). In addition, they can also 
enable a better calculation of the disk thermal physics, 
because the main source of heating for flared disks at 
distances where fragmentation takes place is mostly stel¬ 
lar irradiation (viscosity is more important in the inner 
several tens of AU). 

More details on the coupling of the two numerical codes 
(the hy drodynamic and stellar eyolution ones) can be 
found in lVorobypy et al.l l|2013D . 

2.3. Initial conditions 

For the initial distribution of the gas surface density 
£ and angu lar yelocity fl, we adopted those deriyed by 
iBasd (jl997l l for pre-stellar cores formed as a result of am- 
bipolar diffusion, with the angular momentum remaining 
constant during axially-symmetric core compression 



Here, Hq and £o are the angular yelocity and gas surface 
density at the center of the core and tq = '/Ac^/ttGTio 
is the radius of the central plateau, where Cs is the initial 
sound speed in the core. The gas surface density distribu¬ 
tion described by equation (jO]) can be obtained (to within 
a factor of unity) by integrating the three-dimensional 
gas density distribution characteristic of Bonnor-Ebert 
sp heres with a positiye density-perturbation amplitude 
A (jPapD fc Basd[200^ . In all models the yalue of the 
initial density enhancement A is set to 1.2 and all cores 
haye a fixed ratio rout/?'o = 6.0, where rout is the radius 
of the core, implying that the cores are initially unstable 
to grayitational collapse. 

® Stellar radius affects the proper calculation of the accretion 
luminosity. 


Indiyidual models are generated by first choosing rout 
and then calculating Eg. The central angular yelocity 
Hq is chosen so as to generate cores with different ratios 
of rotational to grayitational energy /3, co nsistent with 
the ya lues inferred for pre-stellar cores by ICaselli et al.l 
(l2002t) . We haye considered 7 models, the parameters of 
which are listed in Table [T] The prototype model 1 has 
the core mass Mcore set to 1.1 Mq, the yalue of /3 set 
to 0.88%, and the yiscous a-parameter set to 5 x 10“®. 
The initial gas temperature is TiuR = 10 K and mag¬ 
netic fields are turned off. Other models are compared 
against model 1, with their parameters being yaried so 
as to emphasize the effect of different initial core masses, 
temperatures, and rotation rates, as well as the effect of 
turbulent yiscosity and magnetic field. 

3. VARIABLE ACCRETION WITH EPISODIC BURSTS 

We start with comparing the long-term eyolution of 
circumstellar disks in our non-magnetic models, continue 
with analysing the accretion rates and the burst statistics 
in different models, and finish with considering the effect 
of magnetic fields. 

3.1. Disk evolution 

Figured] presents the time eyolution of the gas surface 
density £ in seyen models, each row of images corre¬ 
sponding to a particular model. Indiyidual images show 
the inner region of 2000 AU on each side, whereas the to¬ 
tal computational domain is usually much greater. The 
minimum yalue of the density scale is set to -1.5 (in log 
units of g cm~^ ), a typical yalue at the disk outer edge 
(jVorobyoyll2010t) Four columns show the disk eyolution 
at four representatiye times t elapsed since the formation 
of the central protostar. 

We start by analyzing the prototype model 1. Ey- 
idently, the early disk eyolution (t=0.15-0.25 Myr) in 
this model is characterized by yigourous grayitational 
instability and fragmentation. The disk has an irregu¬ 
lar spiral structure with multiple fragments forming in 
the densest parts of the arms. As time progresses to 
t =0.5-0.75 Myr, the disk becomes less irregular, tak¬ 
ing a smoother shape and exhibiting only a week spiral 
structure. No fragments are yisible in the disk at that 
time, which does not necessarily mean that disk fragmen¬ 
tation has ceased completely. As will be shown below, 
this is merely an artifact of infrequent time sampling in 
Figure[I]and short-liyed fragmentation episodes continue 
to occur eyen in the late disk eyolution. 

Neyertheless, the presence of fragments in the early 
disk eyolution and the lack of them in the late eyolu¬ 
tion implies the existence of efficient mechanisms lead¬ 
ing to their loss and/or destruction. Among such 
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Fig. 1.— Images of the gas surface density in the inner 2000 x 2000 AU box in models 1-7 (from top to bottom). Each row represents an 
individual model at four characteristic times after the formation of the central star. The star is marked by the red circle in the coordinate 
center. The scale bar is in log g cm“^ 
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mechanisms are inward migration of fragments onto 
the protostar caused by the gravitational exchange of 
angular momentum with spiral arms a nd other frag¬ 
ments (fVorobvov fc Basul 120061 l2010at iMachida et alJ 
I 2 OIII: iTsukamoto et al.ll2014ll dispersal of fragment s due 
to tidal torques ( Bolev et ^1201(11: IZhu et ani2012ll , and 
ejection of fragments from the disk into the intracluster 
medium due to multi-fra gment gravitational interaction 
(|Basu &: Vorobvo^l2012h . In our case, all three mecha¬ 
nisms were found to be at work. In particular, the inward 
migration and infall of fragments onto the star triggers 
mass accretion bursts discussed in more detail in Sec¬ 
tion 13.21 while the ejection of a fragment with mass on 
the order of O.I Mq at t ~ 0.25 Myr results in a notable 
drop in the net disk mass and a consequent decrease in 
the disk fragmentation activity. 

The Mcore=I-5 Mq model 2 is characterized by a 
greater initial core mass than that of model I. Nev¬ 
ertheless, the evolution in both models is qualitatively 
similar - the disk is vigorously unstable to fragmenta¬ 
tion in the early evolution showing multiple fragments 
interconnected with dense spiral filaments, but becomes 
considerably smoother after t = 0.5 Myr. This transfor¬ 
mation is caused by the same effect in both models - a 
fragment is ejected form the disk leading to an apprecia¬ 
ble drop in the disk mass and subsequent disk stabiliza¬ 
tion. The disk size and mass in model 2 at later times 
are somewhat smaller than those in model 1, but this is a 
mere consequence of a somewhat more massive fragment 
being ejected from the disk in model 2 (~ 0.25 Mq). 

On the other hand, model 3 with a smaller initial core 
mass Mcore=0.3 Mq demonstrates a drastically differ¬ 
ent evolutionary pattern from that of models I and 2 
- the disk exhibits a flocculent spiral structure in the 
early evolution (t < O.I Myr) and becomes increasingly 
axisymmetric with time. No fragmentation is evident 
in the disk (but see Figure [3] below). Model 3 owes 
its special behaviour to the fact that cores with lower 
mass (but with similar /3) form lower mass disks. In the 
case of model 3, the maximum disk mass is 0.09 Mq 
at t Ki O.I Myr and it gradually dr ops to 0.04 Mm at 
t = 1.0 Myr. According to figure 1 in IVorobvovI (l2013ll . a 
model with Mcoie = 0.31 Mq, /3=0.88%, and the disk 
mass of 0.09 Mq lies very near to the fragmentation 
boundary in the /3 — Mcore phase space. This means 
that model 3 can at best produce a couple of isolated 
fragmentation episodes, which might have been missed 
in Figure [T] due to infrequent time sampling. 

Model 4 is designed to highlight the disk evolution re¬ 
sulting from the collapse of a cloud core with low angular 
momentum. In this model, the ratio of rotational to grav¬ 
itational energy /3 is set to 0.14%, more than six times 
smaller than that in model 1, but both models have the 
same value of the initial core mass, Mcore=l-l Mq. As a 
direct consequence of the low initial angular momentum 
in the core, the disk radius in model 4 is rather small 
(ft! 50 AU at t < 0.1 Myr and < 150 AU in the later evo¬ 
lution), much smaller than in other models. No wonder 
that there was only one fragment formed at t = 0.36 Myr 
and that one dispersed after a few orbital period (with¬ 
out reaching the sink cell), possibly due to tidal torques 
from spiral arms. 

Model 5 is set to imitate the effect of a warmer star 
formation environment with a background radiation tem- 


0.0 0.1 


Time (Myr) 



0.3 0.4 

Time (Myr) 

Fig. 2.— Number of fragments vs. time in models 1, 2, 3, 5, 6, 
and 7 (from top to bottom). The number of fragments at a given 
time instant is calculated u s ing th e fragment tracking algorithm 
described in I Vorobyov et ^ ||2()13I1 . An increase in the number of 
fragments shows recent fragmentation, and a decrease shows recent 
destruction/accretion of the fragments. The time is counted since 
the formation of the protostar. 


perature of Tbg = 25 K, 2.5 times higher than that in 
model 1. An increase in Tbg has notably stabilized the 
disk against fragmentation. There are still occasional 
fragmentation episodes taking place in model 5 (see Fig¬ 
ure [2] below), but the disk quickly stabilizes even against 
gravitational instability (put aside fragmentation!) and 
becomes virtually axisymmetric after t = 0.5 Myr. 

The last but one row in Figure [T] presents the disk evo¬ 
lution in model 6 characterized by the a-parameter equal 
to 0.03, six times greater than in model I. The other pa¬ 
rameters are identical in both model 6 and 1. A higher 
efficiency of viscous mass and angular momentum trans¬ 
port does not suppress fragmentation, even though the 
disk mass decreases from a maximum value of 0.25 Mq in 
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model 1 to 0.15 Mq in model 6. Most curiously, model 6 
demonstrates a survival of several fragments orbiting the 
host star at wide separation orbits®. T his is a rare event, 
accordi ng to numerical simul ations of IVorobvov &: Basul 
(|2010bll and IVorobyp-^ (|201.1[ 1 taking place in one out of 
ten models with similar characteristics. We do not think 
that the fragment survival is caused by a higher value 
of a, since the aforementioned studies adopted a smaller 
value of a = 5 X 10“®. Instead, the present simula¬ 
tions suggest that the fragment survival can take place 
for a wide range of a, meaning that this is a robust phe¬ 
nomenon. 

One possible reason why model 6 demonstrated the 
fragment survival, while models 1 and 2 did not, is that 
the latter models experienced fragment ejection, losing 
0.1 Mq and 0.25 Mq of the disk mass, respectively. This 
led to significant weakening of gravitational instability 
and virtual termination of disk fragmentation in the later 
evolution. On the other hand, model 6 lost only 0.02 Mq 
via ejection, which did not affect appreciably the strength 
of gravitational instability and fragmentation. The disk 
in model 6 continued to experience fragmentation even 
after 0.5 Myr of evolution, which greatly increased the 
odds for fragment survival. It appears that the gravi¬ 
tational interaction between fragments in the disk is in¬ 
trinsically a chaotic process, leading in some models to 
fragment ejection and in others to fragment survival. 

Finally, the bottom row in Figure [T] presents the disk 
evolution in model 7 characterized by the non-zero mass- 
to-flux ratio /tb = 3.33. Other parameters in model 7 are 
identical to those of the non-magnetized model 1. The vi¬ 
sual comparison of models 1 and 7 reveals that the frozen- 
in magnetic field does not significantly change the disk 
propensity to fragment: the fragments are present in the 
disk during its early evolution. In the late evolution (af¬ 
ter t = 0.5 Myr), the disk in the magnetized model seems 
to be more extended than its non-magnetized counter¬ 
part, but both show no signs of fragmentation. More 
accurate numerical simulations with non-ideal magneto- 
hydrodynamical effects, such as ambipolar diffusion and 
magnetic braking, are planned for the near future. 

We now analyze the efficiency of disk fragmentation 
in each model using a much higher time sampling than 
in Figure [T] Since we do not introduce sink particles to 
replace fragments in the disk in our Eulerian numerical 
code, it is very difficult to track the position and the fate 
of every fragment during the simulations. Therefore, we 
us ed the fragment detect ion algorithm described in detail 
in IVorobvov et al.l (120131) to postprocess our results and 
calculate the number of fragments present in the disk at 
a given time. We discard fragments that are resolved by 
less than 10 grid cells (3 cells in each direction) since their 
identification on the numerical grid may be dubious. Fig¬ 
ure [2] shows the number of fragments A^f calculated every 
2000 yr after the formation of the protostar in models 1-7 
(from top to bottom). Evidently, Nf varies significantly 
with time and from model to model. An increase in the 
number of fragments indicates recent fragmentation, and 
a decrease implies recent destruction/accretion/ejection 
of the fragments. 

The Mcore = 1.1 Mq model I demonstrates a strong 

® We extended the run time to 1.0 Myr and confirm that the 
fragments were still present in the disk. 


disk fragmentation activity in the early evolution. The 
mean number of fragments in a time period between 
0.03 Myr and 0.26 Myr is four. Episodically, the num¬ 
ber of fragments may exceed 10 or drop to just a few. 
This means that the fragment formation and destruc¬ 
tion mechanisms are constantly at play in the disk. At 
t Ri 0.25 Myr a massive fragment with some circumfrag- 
ment material is ejected from the disk, reducing the to¬ 
tal disk mass and weakening the gravitational instabil¬ 
ity. Part of the ejected material later falls back onto 
the disk triggering two isolated episodes of disk fragmen¬ 
tation at t ~ 0.4 Myr and t r; 0.7 Myr but none of 
those fragments survive. The fragmentation activity in 
the Mcore = 1.5 Mq model 2 is rather similar to that 
in model I except that there are no late-time fragmenta¬ 
tion episodes in the disk, probably due to the fact that 
too much disk mass was lost during the ejection episode 
at t R! 0.32 Myr and little fell back onto the disk. The 
Mcore = 0.3 Mq model 3 shows just a few isolated disk 
fragmentation episodes with the number of fragments 
hardly exceeding one at a time. We do not show model 4 
as it demonstrated no disk fragmentation. 

The disk fragmentation activity in the Tbg=25 K 
model 5 is confined only to the initial O.I Myr of disk 
evolution. The likely explanation is that disk fragmen¬ 
tation in this model is driven by mass infall from the 
collapsing envelope, which is higher than in other mod¬ 
els (Minfaii oc Cg). As Eigure [6] in Section [T2] demon¬ 
strates, Miufaii is maximal at the time of the protostar 
formation and quickly drops afterwards. The disk frag¬ 
mentation process in the q;=0.03 model 6 continues for 
the whole duration of the simulation. We saw in Eig¬ 
ure [T] that several fragments managed to survive through 
almost 1.0 Myr of disk evolution. Figure [3] reveals that 
the number of fragments in the late evolution varies be¬ 
tween one and just a few, meaning that actually only 
one fragment has evolved into a stable companion on a 
wide orbit. Other fragments form in the disk around the 
companion and in the spiral density waves connecting 
the companion with the primary disk rather than in the 
disk around the central star. This interesting effect will 
be studied in more detail in a followup paper. 

Finally, the bottom panel in Figure[2]presents the num¬ 
ber of bursts vs. time in the magnetized model 7. Ev¬ 
idently, the frozen-in magnetic field with a mass-to-flux 
ratio of 3.33 does not appreciably change the disk frag¬ 
mentation activity. The maximum number of fragments 
present in the disk at a specific time = 8) is some¬ 

what smaller in model 7 than in the corresponding non- 
magnetized model I = 12), but otherwise the time 

behaviour of iVn is similar in both models. 

3.2. Accretion and infall rates 

In this section, we analyze the time behavior of mass 
accretion rates onto the star (M) and infall rates onto 
the disk (Minfaii) obtained in models 1-6. We note that 
M is calculated at the position of the inner sink cell, 
Tcs = 6 AU, while Minfaii - at a distance of 2000 AU. 

FigureOpresents M (black solid lines) and Minfaii (red 
dashed lines) vs. time in models 1-3, which are char¬ 
acterized by different initial core masses Mcore as indi¬ 
cated in each panel. Other parameters in these mod¬ 
els are similar. The time is counted from the begin- 
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ning of numerical simulations (i.e., from the beginning 
of the core contraction). The formation of the cen¬ 
tral protostar in models 1-3 occurs around 0.15 Myr, 
0.2 Myr, and 0.05 Myr, respectively. Evidently, accre¬ 
tion variability increases with increasing Mcore- While 
the Mcore = 0.3 Mq model 3 shows only an order of 
magnitude variations in M with just a few stronger fluc¬ 
tuations, the Mcore = 1-1 Mq and Mcore = 1-5 Mq 
models 1 and 2 are characterized by highly variable ac¬ 
cretion spanning several orders of magnitude with peak 
values exceeding in magnitude 10“^ Mq yr“^. These 
accretion bursts are caused by fragments spiralling into 
the central star due to the loss of angular momentum 
via gravitational interact ion with spiral arms and other 
fragme nts in the disk (|Vorobvov fc Basul 120051 120061 
[20103)- The animation of this process can be viewed at 
http : //www. astro .uwo . ca/^basu/movies .html. We 
note that in model 2 one accretion burst at ~ 0.52 Myr is 
especially strong, exceeding in magnitude 10“^ Mq yr“^. 
This burst is coeval with the fragment ejection event dis¬ 
cussed in the previous section, a pair phenomenon that 
causes one fragment to fly into the intracluster medium 
and the other to fall o nto the star (for more detail see 
iBasu fc Vorobvo\^l2012ll . 

We now discuss the reason why the time behaviour of 
mass accretion rates in models with increasing Mcore is 
so different. Figure [T] has already given us a hint - gravi¬ 
tational instability in models 1 and 2 is notably stronger 
than in model 3. This is not surprising considering that 
more massive cores are supposed to form more massive 
disks. To quantify the strength of gravitational insta¬ 
bility in the disks of our models, we calculated global 
Fourier amplitudes using the following equation: 


Crrr{t) 


1 



cRd 


E(r, (/), t) dr d(j) 


( 11 ) 


where M^ is the disk mass, i?d is the disk’s physical outer 
radius, and m is the number of the spiral mode. When 
the disk surface density is axisymmetric, the amplitudes 
of all modes are equal to zero. When, say, Cm{t) = 0.1, 
the perturbation amplitude of spiral density waves in the 
disk is 10% that of the underlying axisymmetric density 
distribution. 

The Fourier amplitudes of the first three spiral modes 
are presented in Figure H] for models 1-3. Evidently, the 
Mcore=0.3 Mq model 3 is characterized by the lowest 
Fourier amplitudes - they hardly exceed 10% that of the 
underlying axisymmetric density distribution in the early 
evolution and quickly decline with time. The lower-order 
modes are generally higher in amplitude than the higher- 
order modes. On the other hand, the amplitudes of spiral 
modes in models 1 and 2 are appreciably higher than in 
model 3, reaching values as high as 60% that of the un¬ 
derlying axisymmetric distribution (log Cm ~ —0.2) by 
the end of the embedded phase. In the later evolution, 
the Fourier amplitudes gradually decline with time, re¬ 
flecting the overall disk stabilization due to continuing 
loss of disk material via accretion onto the star. 

When comparing Figures [3] and 01 we can notice a gen¬ 
eral correlation between Fourier amplitudes and variabil¬ 
ity in the mass accretion rates. The low-M^ore model 3 
is characterized by both low Fourier amplitudes and 
low accretion variability. As Cm declines with time, 
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Fig. 3.— Mass accretion rates at 6 AU (black solid lines) and 
envelope infall rates at 2000 AU (red dashed lines) in models 1-3. 
The arrows mark the formation of the first hydrostatic core and 
the disk in model 1. 

the variability in M diminishes. For Cm 0.01, the 
non-axisymmetric density perturbation in the form of 
spiral density waves is only 1% that of the underly¬ 
ing axisymmetric distribution, implying that the mass 
transport is now mostly controlled by viscous torques 
(jVorobvov fc Basi] l2009h . The latter drive the disk 
toward an axisymmetric state as evident in Figure [T] 
and the corresponding accretion rates show only low- 
amplitude flickering. On the other hand, Fourier am¬ 
plitudes in the higher-Mc models 1 and 2 are greater 
than in model 3 and the variability in the correspond¬ 
ing accretion rates is notably higher. Fourier amplitudes 
comparable to or greater than 10% imply the presence 
of strong non-axysimmetry in the disk in the form of 
both spiral waves and fragments. It is important to note 
that Fourier amplitudes are highly variable^ due to the 
dynamical interaction of spiral density waves and frag¬ 
ments in the disk. A combination of two effects: that of 

^ In fact, the time variability in Cm is even higher but we had 
to smooth Fourier amplitudes somewhat to make them discernible 
in the figure. 
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Fig. 4.— Global Fourier amplitudes of the first three spiral modes 
m=l-3 in models 1-3. The horizontal dotted lines mark the bound¬ 
ary above which spiral waves have perturbation amplitudes greater 
than 10% that of the underlying axisymmetric density distribution. 

the nonlinear interaction between different spiral modes 
and fragments in the disk and that of the fragments spi¬ 
ralling down onto the star produces variable accretion 
with episodic bursts. 

We now proceed with analyzing the effect of different 
values of j3, Tbg, a-parameter on the variability of M. 
Figure [5] presents the mass accretion and infall rates in 
model 4 characterized by /3 = 0.14%, a factor of six lower 
than in model 1. The other parameters in both models 
are identical. Evidently, the mass accretion rate in the 
low-/? model 4 is characterized by low-amplitude flick¬ 
ering and complete absence of bursts, while the high-/? 
model 1 shows strong accretion variability with multiple 
bursts, some exceedin g in magn i tude 1 0~^ Mq yr“^. As 
was demonstrated in IVorobvovI (|2010ll . pre-stellar cores 
with higher /? tend to form more massive and extended 
disks than cores with lower /?, which can be understood 
on the basis of simple centrifugal radius arguments. Mas¬ 
sive and extended disks are more gravitationally unstable 
and prone to fragmentation than light and compact ones, 
explaining the aforementioned tendency in M. 

As a next step, we describe the effect of a higher back- 



Fig. 5.— Mass accretion rates at 6 AU (black solid lines) and 
envelope infall rates at 2000 AU (red dashed lines) in models 4. 


ground temperature Tbg on the time behaviour of mass 
accretion rates, thus mimicking a higher heating rate 
coming from the external environment. The top panel 
in Figure [6] presents M and Minfaii vs. time in model 5 
characterized by Tbg = 25 K, which is 2.5 times higher 
than the corresponding value in model 1. The other pa¬ 
rameters are identical in both models. The increase in 
Tbg reduces notably the accretion variability. This effect 
is explained by the fact that a higher background tem¬ 
perature raises the overall disk temperature and reduces 
the strength of disk gravitational instability. Neverthe¬ 
less, the Tbg=25 K model 5 exhibits several bursts with 
M ^ 10-4 Mq yr-i. 

The middle panel in Figure [6] presents M and Minfaii 
vs. time in model 6, the parameters of which are similar 
to those of model 1 except that the a-parameter is set 
now to 0.03, a factor of six higher than in model 1. Ev¬ 
idently, the increase in a acts to reduce the variability 
in M, which now features only order-of-inagnitude vari¬ 
ations. Nevertheless, there are two well-pronounced ac¬ 
cretion bursts with M approaching lO-^ Mq yr”^. The 
notable reduction in the burst frequency in model 6 as 
compared to the low-a model 1 can be attributed to in¬ 
creased viscous mass transport through the disk, which 
reduces both the d isk mass and the strengt h of gravi¬ 
tational instability (|Vorobvov fc Basull2010all . Viscosity 
also tends to s mooth out local non-axisy mmetric density 
enhancements (|Vorobvov fc Basi] 120091 ) . thus reducing 
the accretion variability caused by spiral density waves. 

The last model in this study takes into account the 
effect of frozen-in magnetic fields. The bottom panel 
in Figure [5] presents the mass accretion and infall rates 
in model 7 which has the same parameters as the non¬ 
magnetic model 1 except that the mass-to-flux ratio is 
set to = 3.33. Evidently, the magnetized model 7 
exhibits several strong accretion bursts with a rate > 
10-4 Mq yr-4, indicating that the frozen-in magnetic 
field does not suppress the burst phenomenon. 

Finally, we briefly discuss the time behaviour of in¬ 
fall rates in our models and a possible link between 
Afinfaii and the burst phenomenon. Figures [H [5l and 
[6] demonstrate that Minfaii steadily increases from « 
10-® Mq yr-4 to (5 - 10) x IQ-® Mq yr-^ during 
the early evolution and then gradually declines below 
10-’’ Mq yr-4, reflecting the overall depletion of mass 
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Fig. 6.— Mass accretion rates at 6 AU (M, black solid lines) 
and envelope infall rates at 2000 AU (Minfaii, red dashed lines) in 
models 5, 6, and 7 (from top to bottom). 


in the parental cores. Order-of-magnitude variations of 
Ailinfaii in the late evolution of models 1, 2, and 6 are 
caused by the disturbing influence of fragments that were 
ejected from the disk. A visual inspection of M and 
Vfinfaii in Figures [3l [3 and |6] reveals that the amplitude 
of variations in the accretion rate generally correlates 
with the infall rate, though with some notable excep¬ 
tions. More specifically, for 10“® Mq yr“^ < Minfaii < 
10“® Mq yr“^ both the large-scale accretion variabil¬ 
ity and the bursts are present. For 10“^ Mq yr“^ < 
Afinfaii < 10“® Mq yr“^ some variability in accretion 
rates is still present but the bursts are mostly gone, ex¬ 
cept for model 6 showing one energetic burst in the late 
evolution. For Minfaii < 10“^ Mq yr“^ the accretion 
variability diminishes and no bursts are seen. The char¬ 
acter of accretion for various infall rates is summarized 
in Table [3 

Model 4, however, stands apart and shows little accre¬ 
tion variability and no bursts even though the infall rates 
are greater than 10“® Mq yr“^ during the early evolu- 


TABLE 2 

Infall rates and the character of accretion 


^infali Variability Bursts 

(Mq yr“^) (orders of mag.) 

10-®-10-“ 2^^3 multiple 

10“®-10“^ 1 occasional 

< 10“^ <1 no 


Time (Myr) 

0.250 0.255 0.260 0.265 0.270 0.275 



Fig. 7.— Comparison of the mass accretion rates at 6 AU and 
8 AU (top panel) and 6 AU and 12 AU (bottom panel) during a 
short period of evolution in model 1. 

tion. Due to its low angular momentum, the /3 = 0.14% 
model 4 starts forming the disk considerable later than 
other models. As a result, the disk mass and radius are 
not sufficiently large for gravitational fragmentation to 
take place (see also Fig[T|). High infall rates are therefore 
a necessary but not sufficient condition for the develop¬ 
ment of the burst phenomenon. 

To summarize this section, our numerical simulations 
with the updated code con fi rmed t he findings reported 
earlier in IVorobvov fc Basul (|2010aD , namely that an in¬ 
crease in the initial core mass M^ote and/or the ratio 
of rotational to gravitational energy j3 acts to increase 
the amplitude of accretion variability and the number of 
accretion bursts. An increase in the a-parameter and 
Tbg is found to have the opposite effect - the accretion 
variability and bursts diminish but do not cease to ex¬ 
ist, at least for reasonable values of a and Tbg. This 
means that accretion variability and bursts are a robust 
phenomenon, weakly sensitive to (modest) variations in, 
e.g., dust opacities, turbulent viscosity, and stellar radi¬ 
ation. 

3.3. The effect of the inner boundary condition 

In our models, the mass accretion rate M is calculated 
at the position of the inner sink cell, r^c = 6 AU. Two 
questions arise in this context: how much M is sensitive 
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both to the choice of Tsc and to the imposed free outflow 
boundary condition through the sink cell. We cannot 
compute the disk dynamics at sub-AU scales because of 
a strict Courant time step limitation imposed on our ex¬ 
plicit Eulerian code. However, we performed several test 
runs with rgc varied by a factor of 2 and found little effect 
on the accretion burst phenomenon. 

This result can be easily understood by the following 
reason. Fragments form in the disk at a radial distance 
of at least several tens AU from the sink cell where condi¬ 
tions beconie_^vourable for gravitational fragmentation 
(e.g. iMeru fc Batel[2?)T^ . At such distances, the influ¬ 
ence of the inner boundary on the disk fragmentation 
process is negligible. Furthermore, the typical size of 
fragments approaching the inner sink cell is compara¬ 
ble to the size of the sink cell itself. Therefore, vary¬ 
ing Tsc by a factor of 2 does not make much difference 
- the fragment will sooner or later pass through the 
sink cell. On the other hand, as fragments approach 
the star, they must be inevitably stretched out due to 
tidal torques. How much of the fragment material fi¬ 
nally reaches the star is an op en question and requires 
a focused investigation (se e e.g. ICha fc Navakshinll2nill : 
iNavakshin fc Lodat3l2012D . The effect of a sudden mass 
deposition onto the inner disk (< 10 AU), as if by infall 
of a fragment migrating through the disk onto the s tar, 
has recently been investigated bv lOhtani et al.l l)2014[) . It 
was found that such an event can lead to the FU-Orionis- 
like eruption due to triggering of the magneto-rotational 
instability (MRI) at sub-AU scales. This means that ei¬ 
ther directly, by deposition of material onto the star, or 
indirectly, by triggering the MRI in the inner disk, the 
inward migration of fragments will likely produce lumi¬ 
nosity outbursts. 

To evaluate the possible influence of the free outflow 
boundary condition on the mass accretion rate through 
the sink cell, we calculated the mass transport rates at 
a few AU away from the sink cell. Figure [7] presents the 
mass accretion rates through the sink cell at 6 AU (M, 
black lines), as well as the mass transport rates through 
the disk at 8 AU {Ms, red line) and 12 AU {M 12 , red 
lines). A narrow time interval of 2 x 10^ yr in model I 
is chosen to focus on a time period featuring both the 
bursts and quiescent accretion. Evidently, Mg is very 
similar to M except for a few instances when notable 
deviations are visible. On the other hand, M 12 demon¬ 
strates larger deviations from M, but retains the main 
qualitative features of M such as accretion bursts. Both 
Mg and Mi 2 exhibit more variability than M, most likely 
due to the fact that the inner boundary allows for matter 
to flow into the sink cell but not out of it, thus somewhat 
artificially damping the time variations. 

4. CHARACTERISTICS OF LUMINOSITY BURSTS 

In this section, we analyze the characteristics of lu¬ 
minosity bursts obtained in our models and compared 
them with the available statistics on FU-Orionis-type 
eru ptions (FUors ) taken from the recent review paper 
by lAudard et al.l (|20I4ll . In order to distinguish the 
bursts from regular (order-of-magnitude) variability in 
our models we have to make several assumptions. First, 
we assume that the total luminosity L* during the burst 
should be comparable to that of FUors. The latter are 


usually characterized by an increase in brightness by at 
least a factor of 3-4 (in stellar magnitudes) as compared 
to the pre-burst, quiescent phase. Therefore, we stipulate 
that the luminosity increase during the burst should be 
at least 16 times 3 mag) that of the pre-burst phase. 

Calculating the luminosity in the pre-burst phase 
turned out to be not an easy task due to a highly change¬ 
able nature of accretion. We do this by defining the so- 
called background luminosity Tbg, which comprises the 
stellar photospheric luminosity and the mean ac¬ 

cretion luminosity (T*,accr)- The former is provided by a 
stellar evolution code (see Section [2|), while the latter is 
found as 


(T^^accr) — 


GM*(M) 
2i?* ^ 


( 12 ) 


where (M) is the mean accretion rate calculated using a 
running average of the instantaneous accretion rates M 
over a time period of 10^ yr. When doing the average, we 
filtered out values that are greater than 5 x I0“® Mq yr“^ 
by the reason that they may already represent a burst in 
its rising or fading phase. 

The red and black lines in Figure |8] present the total 
luminosity A* and background luminosity Abg in mod¬ 
els I, 2, 5, and 7. We left out models that showed too 
few bursts to be statistically meaningful. In general, A* 
is highly variable in the early evolution, reflecting the 
corresponding variations in the mass accretion rate. On 
the other hand, Abg shows much less variability and de¬ 
scribes well the minimal luminosity in each model. The 
blue lines mark the values that are 16 times greater than 
the background luminosity at a given time, representing 
therefore the 3-magnitude cutoff above which a surge in 
luminosity may be classified as a FUor. In principle, 
these relatively modest bursts can be confused with the 
so-c alled EXors named after its prototype EX Lupi (see 
e.g. lAudard et al.ll20I4ll . Therefore, with the blue line 
we also plot the 4-magnitude cutoff (39 times greater 
than Abg) in order to analyze the statistics of more ener¬ 
getic bursts, which are more likely to represent bona fide 
FUors. 

Evidently, models I and 2 are characterized by the 
largest number of strong FUor-type bursts, amounting 
to 10 and more per model®. The higher-Abg model 5 has 
only a few bursts above the 4-magnitude cutoff. At the 
same time, model 7 demonstrates several strong bursts 
despite the presence of frozen-in magnetic field with a 
mass-to-flux ratio /ie = 3.33, indicating that FUors can 
occur in magnetized disks as well. 

In the following text, we analyze the main character¬ 
istics of the luminosity bursts obtained in our models 
in order to compare our predictions with observations. 
Figure |9] presents the duration of the burst fbst (left col¬ 
umn), the accreted mass during the burst Maccr (middle 
column), and the peak luminosity during the burst A((^“ 
(right column) in models I, 2, 5 and 7 (from top row to 
bottom one). Only bursts with a 4-mag cutoff are shown. 
The x-axis shows the ordinary number of the burst ar¬ 
ranged along the line of the burst occurrence. Evidently, 
the burst duration stays mostly in the 10-100 yr limit 
with little dependence either on time or particular model. 


® Some of the bursts are closely packed and cannot be resolved 
in the figure (see Section (6)1. 
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Fig. 8.— Red lines: total (accretion plus photospheric) luminosity vs. time in model 1 (top left), model 2 (bottom-left), model 5 
(top-right), and model 7 (bottom-right). The black lines provide the background luminosity comprising the photospheric luminosity plus 
accretion luminosity arising from accretion with a rate < 5 X 10“® Mq yr“^. The blue and pink lines mark the 3-magnitude and 4- 
magnitude cutoffs above which a surge in luminosity is considered to be an FU-Orionis-type outburst. The vertical dotted lines mark the 
Class O/I boundary (left lines) and Class I/II boundary (right lines). See the text for more details. 


These values are in good agreement with the measured 
or inferr ed duration of FUor s (see Table [2] below and ta¬ 
ble 1 in lAudard et al.l[2014D . The accreted mass during 
the bursts ranges from 1.0 to 75 Jupiter masses, covering 
the full mass range of giant planets and brown dwarfs 
and refle cting the mass range of fragments forming in 
the disk (|Vorobvov et al.ll2013l) . The peak luminosities 
of most bursts span a range between 75 Lq and 600 T©, 
with a few notable exceptions in model 2 reaching values 
in excess of 3000 L©. These most luminous (and closely- 
packed) bursts occur during the fragment ejection event 
discussed in Section o and are the result of the conser¬ 
vation of angular momentum, causing one fragment to fly 
out of the disk and the other fragment to fall onto the star 
due catastrophic loss of angular momentum during the 
close encounter. The fact that there are three closely- 
packed bursts instead of just one is explained by tidal 
destruction of the infalling fragment (see Section [H). 

The summary of various characteristics of the bursts 
obtained in our modeling are provided in Table [31 More 
specifically, the second column provides the total num¬ 
ber of bursts iVbst and the number of bursts in the deeply 
embedded Class 0 phase (in parentheses, see Section [5]), 
the third column is the fraction of stellar mass accreted 
during the bursts and the fourth column is the 

fraction of total disk lifetime spent in the burst phase 
^bst- other columns (from 5th to 9th) present the 

maximum, minimum and mean luminosities of the bursts 
{Lmax, Amin and Amean), the maximum, minimum and 
mean accretion rates during the bursts 
and Mmean), the maximum, minimum and mean burst 
durations the maximum, mini- 



1 2 3 

Burst number 


12 3 12 3 

Burst number Burst number 


Fig. 9.— Burst characteristics in model 1 (top row), 2 (middle 
row) and 5 (top row). Columns (from left to right) present the 
duration of the bursts tbst (ii^ y^)? the accreted mass during the 
bursts Maccr (in Mj^p), and the peak luminosity during the bursts 
(” Lq). 


mum and mean duration of the quiescent phase between 
the bursts® and and the maximum 

® We calculated these quantities by making no distinction be- 
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and mean accreted mass during the bursts and 

The mean values were found by arithmetically 
averaging over all bursts in a particular model. The 
known characteristics of FUors are provid ed in the bot¬ 
tom l ine and are taken from table 1 of lAudard et al.l 

(IMl . 

Our models seem to reproduce the main properties 
of FUors rather well. When averaged over all models, 
the mean luminosity during the burst in our models is 
312 Lq, about a factor of 1.5 greater than that of FUors. 
The agreement becomes even better after taking out rare 
outliers such as very energetic bursts in model 2. The av¬ 
eraged burst duration is 39 yr, again almost a factor of 2 
greater than that of FUors, but this may be simply due 
to the fact that many FUors are still found in the active 
phase. The duration of the quiescent phase also seems to 
be in agreement with the recently estimate d lower limit 
of 5-1 0 kyr using mid-infrared photometry (jScholz et al.l 
l2013ll . There is however some disagreement which may 
be of physical nature. For instance, our model accretion 
rates agree well with the mean inferred mass accretion 
rates in FUors, but fail to explain low-M objects such as 
HBC 722 with M ~ 10“® Mq yr“^ (lAudard et al.l[20l^ . 
Indeed, for the minimum accreted mass of 1.0 Jupiter and 
the maximum duration of the burst of 100 yr, the result¬ 
ing minimum accretion rates in our models are supposed 
to be around 10“® Mq yr“^, implying that FUors like 
HBC 722 may be driven by mechanisms other than disk 
fragmentation. 

5. EMBEDDED VS. OPTICALLY VISIBLE BURSTS 

According to iQuanz et al.l (|2007ll . FUors can be clas¬ 
sified in two categories, depending on whether silicate 
features at 10 /rm are seen in absorption or emission. 
FUors with silicate in absorption are likely still embed¬ 
ded in parental envelopes, whereas FUors with silicates in 
emission are likely more evolved and having (partially) 
depleted envelopes. Out of 21 observed FUors in the 
Quanz et al. sample, 12 have silicates in absorbtion and 7 
in emission, with 2 objects having a flat spectrum (mak¬ 
ing their classification dubious). This simple analysis 
suggests that most FUors are rather young objects, pos¬ 
sessing sizeable envelopes. 

A similar conclusion can be made using data summa¬ 
rized in the r ecent review on episo dic accretion in young 
protostars bv lAudard et al.l (|2014ll . Two prominent ob¬ 
jects, FU Ori itself and V1515, are often considered as 
most evolved FUors having little-to-no envelope, sug¬ 
gested by weak far-infrared/submillimeter continuum be¬ 
yond 100 /im in their spectral energy distributions^'^. Ac¬ 
cording to table 1 in Audard et ah, these two objects have 
the optical extinction Ay lying in the 1.5-3.2 range. On 
the other hand, young FUors with silicate i n absorption 
have a minimum extinction of Ay = 4.2 (| Quanz et al.l 
1200711 . We therefore set Ay = 4.0 as a tentative bound¬ 
ary between embedded and optically visible FUors. In 
the Audard et al.’s sample, 15 FUors are characterized 
by Ay > 4.0 and only 6 FUors have Ay < 4.0. Of 
course, this simple analysis may somewhat be affected 

tween isolated and clustered bursts, see Section [5] 

Even in this case, near-infrared interferometry shows that, 
e.g., V1515 may have an u ndetected contribution from the envelope 
IIMillan-Gabet et ahll^OGfl . 


by inclination: highly inclined FUors may appear em¬ 
bedded, whereas in reality they are not. Nevertheless, 
all above arguments taken together indicate that many 
(if not most) FUors are young, embedded objects rather 
than older, Class II stars. 

To classify the bursts in our models, we use the remain¬ 
ing mass in the envelope to define the boundary between 
the embedded and optically visible phases. Namely, we 
assume that the optically visible Class II begins when 
less than 10% of the initial core mass is left in the enve¬ 
lope. The boundary between the deeply embedded Class 
0 phase and the partly embedded Class I phase is defined 
as the time when 50% of the initial core mass is left in 
the envelope. Our adopted classification scheme is based 
on physical properties of a youn g stellar object, such 
as envelope and disk mass es ('e.er.. iRobitaille et al.ll2006l : 
iDunham et al.l[2010L l20I4ll , rather than on observational 
signatures, such as sub millimeter luminosi t ies or effec¬ 
tive temperatures (e.g. lAndre et al.l 119931 : IChen et all 
[TMl). Classifications relying upon physical properties 
are usually referred in the literature as “stages”, whereas 
those using observational signatures are called “classes”. 
For simplicity here we use the term “class” to refer to 
both the physical stages and observational classes. Our 
adopted defin ition of physica l stage s was extensively in¬ 
vestigated in iDunham et al.l (1201011 . They found that 
there is not always a one-to-one correspondence between 
physical stage defined by the envelope mass and obser¬ 
vational class defined by the submillimeter luminosity or 
effective temperature due to the effects of geometry and 
extinction. In reality the exact point at which to set 
the class boundaries is somewhat uncertain, which could 
shift the duration of the embedded phase in our models 
by a factor of order unity in either direction. We disen¬ 
tangle the disk and infalling envelope on our numerica l 
grid using the algorithm described in IVorobvo~^ (j2011ll , 
which is based on the disk-to-envelope transition density 
of Scrit = 0.5 g cm“^ and the velocity field in the in¬ 
falling envelope. Varying the value of Scrit by a factor 
of 5 results in changes of the estimated onset time of 
different phases by only a few per cent. 

The vertical dotted lines in Figure[8]mark the Class O/I 
boundary (left lines) and the Class I/II boundary (right 
lines). Evidently, most bursts in our models occur in the 
partly-embedded Class I phase. For instance, model I 
has nine strong bursts (above 4-mag cutoff) taking place 
in the Class I phase, and only one in the Class 0 phase, 
while model 2 has only five strong bursts out of 13 oc¬ 
curring in the Class 0 phase. The values in parentheses 
in the second column of Table [3] provide the number of 
bursts (out of the total number) occurring in the Class 0 
phase. Out of the total 29 strong bursts, only 8 occurred 
in the deeply embedded Class 0 phase and 21 in the 
partly-embedded Class I phase. As similar though less 
pronounced tendency is found for less energetic bursts 
(3-mag cutoff), except for the magnetized model 7, in 
which most bursts occur in the early Class 0 phase. 

Our models notably lack bursts taking place in the 
optically visible Class II phase. This is not surprising 
since most fragments form in the early evolution, which 
is characterized by most massive and gravitationally un¬ 
stable disks (see Figs. [T] and [2]). This phase is also less 
favourable for the survival of fragments owing to strong 
gravitational and tidal torques which tend to drive frag- 
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TABLE 3 

Characteristics of luminosity bursts 


Model 

■(Vbst 


.ftot 

^bst 

-^'max/ -^min/ -^mean 

-^max/ '^min/-^mean 

xmax /amin /j.mean 
^bst /^bst /^bst 

j.max /j.min /xmean 
'^qst / ^qst / ‘'qst 

Ti^max /Ti^mean 



(%) 

(%) 

{Le) 

(10-4 Mq yr-1) 

(yr) 

(104 yr) 

(Mjup) 

4-mag cutoff 

1 

10(1) 

3.8 

0.035 

357/87/208 

2.4/0.78/1.4 

48/12/25 

10/1.6/4.7 

4.0/2.3 

2 

13(5) 

18.6 

0.06 

3042/77/846 

20/0.8/5.3 

78/10/36 

15A.36/4.4 

75/12 

5 

3(0) 

1.3 

0.02 

620/403/500 

1.2^.92/1.0 

82/18/41 

- 

4.4/2.5 

7 

3(2) 

2.0 

0.04 

227/115/157 

2.6/0.9/1.5 

64/36/47 

- 

4.1/3.3 

3-mag cutoff 

1 

21(9) 

6.4 

0.086 

357/28/110 

2.4/0.25/0.81 

120/2/31 

4.3/0.005/1.0 

9.0/1.9 

2 

42(15) 

24.7 

0.24 

3042/25/300 

20/0.14/2.1 

320/2/47 

6.4/0.002/0.8 

75/6 

5 

15(0) 

5.2 

0.09 

620/163/290 

1.2/0.21/0.64 

132/14/50 

1.4/0.004/0.33 

5.8/2.4 

7 

12(8) 

4.4 

0.1 

227/37/81 

2.6/0.37/0.79 

90/6/33 

6.2/0.002/1.1 

6.7/1.9 






FUors 









(observations) 





26 

- 

- 

525/10/200 

10/0.01/1.9 

80/4/20 

- 

- 


merits onto the star or destroy them. On the other hand, 
fragments that happen to survive through the embedded 
phase are more likely to form stable companions, rather 
than to migrate onto the star and trigger a burst. 

There are however exceptions. The middle panel in 
Figure [6] shows the mass accretion rate in model 6. The 
Class II phase in this model starts at t = 0.47 Myr after 
the onset of gravitational collapse. A strong accretion 
burst at t = 0.763 Myr, corresponding to a luminosity 
outburst of 373 Lq, occurs well into the optically visible 
Class II phase. 

Model 6 is unique among other models in the sense 
that it reveals the formation of a wide-orbit companion 
on a quasi-stable orbit. Figure [10] presents the gas sur¬ 
face density in model 6 in the inner box of 1400 x 1400 AU 
just before the luminosity outburst and immediately af¬ 
ter it. Evidently, the survived companion possesses a 
circumfragment disk, which is sufficiently massive to ex¬ 
perience episodic fragmentation. Indeed, the mass of the 
companion and its disk are 57 Mjup and 21 Mjup, making 
the disk to central object mass ratio equal to ^ ~ 0.38. 
Systems wi th £ > 0.1 are likely to b e unstable to frag¬ 
mentation ([Vorobyov fc Basul l2010all . For comparison, 
the masses of the central star and its circumstellar disk 
at this time are 0.7 Mq and 0.022 Mq, and the disk 
to star mass ratio is only Ri 0.03, explaining the lack of 
fragmentation in the circumstellar disk. The luminosity 
outburst in model 6 is caused by one of the fragments 
(shown by the yellow arrow) forming in the circumfrag¬ 
ment disk and falling onto the star owing to a complex 
interplay and exchange of angular momentum with other 
fragments and spiral filaments. 

Finally, we note that late bursts in our model are 
also possible if a fragment is ejected from the disk into 
the intrucluster med ium through the multi- body gravi¬ 
tational interaction ([Basil fc Vorobvovll2ni2h . The ejec¬ 
tion is paired with another fragment losing its angu¬ 
lar momentum and falling onto the star, producing a 
strong accretion burst. As table 1 in iBasii fc VorobvovI 
([2012f) demonstrates, some ejection events may occur 
0.7-0.8 Myr after the formation of the protostar, which 
is usually an optically visible phase. To summarize, most 
burst events in our model take place in the partly em¬ 
bedded Class I phase, with a smaller fraction occurring 
in the deeply embedded phase and a few bursts in the 



Fig. 10. — Gas surface density images (in log g cm“^) showing 
the evolution in model 6 before and after the accretion burst at 
t = 0.763 Myr. The yellow arrows track the position of a fragment 
forming in the disk of the companion and migrating into the cen¬ 
tral star. The time is counted since the beginning of gravitational 
collapse. The central star is formed at 0.146 Myr. 


optically visible Class II phase. 

The tendency of luminosity bursts to mainly occur in 
the embedded phase can be understood by analyzing the 
updates applied to our numerical hydrodynamics code. 
The use of Semenov opacities (instead of Bell & Lin’s) 
and stiffer equation of state both contribute to make disk 
fragmentation more difficult owing to an increased disk 
temperature. However, the envelope infall onto the disk 
in the embedded phase increases the surface density and 
brings periodically the disk to the fragmentation bound¬ 
ary. Once the envelope begins to dissipate and the in¬ 
fall rate drops, the burst activity subsides, explaining 
the scarcity of bursts in the Class II phase. The dif¬ 
ference with our previous work here is that it now takes 
pre-stellar cores with somewhat higher mass and angular 
momentum to trigger the burst phenomenon after form- 
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Fig. 11.— Time evolution of the stellar radius (left column), pho- 
tospheric luminosity (middle column) and ratio of total luminosi¬ 
ties T«,Lvon/^*. DAM (right column) deriv ed from the Lyon stellar 
evolution code lIBarafie fc Chabrieil I20 1CII~I and the non- accreting 
stellar evolution tracks of lD’Antona^ Mazzitellil ll 19941V Shown 
are results for model 1 (top row) and model 2 (bottom row). The 
vertical dotted lines mark the Class O/I and Class I/II boundaries. 

ing a star-disk system. 

The use of the Lyon stellar evolution code (instea d 
of precalculated tracks of iD’Antona &: Mazzitellil (|1994f) ) 
has a more complicated effect. The left and middle 
columns in Figure [TT] present the time evolution of the 
stellar radius i?* and photospheric luminosity L*,ph in 
model 1 (top row) and model 2 (bottom row) derived 
using the Lyon stellar evolution code (solid lines) and 
D’Antona & Mazzitelli stellar evolution tracks (dashed 
lines). The right column shows the ratio of total lumi¬ 
nosities L*,Lyon/T*,DAM found using the Lyon code and 
D’Antona & Mazzitelli tracks (hereafter, DAM tracks). 
The vertical dotted lines mark the Class O/I and Class 
I/II boundaries. 

Evidently, the time evolution of A* is different in the 
Lyon code and DAM tracks. In the early Class 0 phase, 
the stellar radius in the Lyon code is significantly smaller 
than that of the DAM tracks. As a result, the accretion 
luminosity in the Lyon tracks is considerably greater^^. 
At the same time, the photospheric luminosity in the 
Lyon code and DAM tracks are similar. The net result 
is that the total luminosity in the Lyon code is notably 
greater in the early Class 0 phase, as the right-hand-side 
column in Figure [TT] demonstrates, making disk fragmen¬ 
tation more difficult and reducing the burst activity in 
the early Class 0 phase. 

In the late Class 0 phase, the stellar radius in the Lyon 
code increases owing to absorption of a fraction of the 
accretion energy and the situation reversers: i?* in the 
Lyon code becomes larger than that in the DAM tracks. 
Steep episodic rises in i?* seen in model 1 and especially 
in model 2 at t « 0.55 Myr are caused by accretion 
bursts, leading to stellar bloating due to the absorbed 
accretion energy. In the Class I phase, A* in the Lyon 
code is systematically larger than in the DAM tracks and 
the photospheric and accretion luminosities are smaller 
(apart from time instances with strong bursts). The net 
result is that the total luminosity in the Lyon code be¬ 
comes smaller on average that that in the DAM tracks, 

Note that when deriving the stellar radii and photospheric 
luminosities from the DAM tracks we use the accretion histories 
obtained by hydrodynamical simulations coupled with the Lyon 
stellar evolution code. That is why the accretion rate and stellar 
mass are identical in both the DAM and Lyon cases. 


Time (Myr) Time (Myr) 



0.2524 0.2526 0.2528 0.2530 0.3888 0.3890 0.3892 0.3894 0.3896 



Time (Myr) Time (Myr) 


Fig. 12.— Zooming in onto individual bursts in model 1 (left 
column) and model 2 (right column). The solid lines are the total 
luminosity vs. time, while the dashed and dash-dotted lines are 
the 4-magnitude and 3-magnitude cutoffs. 


making disk fragmentation easier. This effect, along with 
continuing infall from the envelope, explain why luminos¬ 
ity bursts tend to occur in the Class I phase. 


6. ISOLATED AND CLUSTERED LUMINOSITY BURSTS 

In this section, we zoom in onto individual luminosity 
bursts in our models in order to describe their evolution 
on short time scales. Figure [T^l presents the total lumi¬ 
nosity of individual bursts in model 1 (left column) and 
model 2 (right column) during six time periods, each of 
10^ yr in duration. Also shown with the dashed and 
dash-dotted lines are the 3-magnitude and 4-magnitude 
cutoffs to help identify the bursts (see Section |T|). Two 
different types of bursts are evident in the figure. The 
first type can be described as single, isolated bursts as 
shown in panels c, d, and e. These bursts are caused by 
compact infalling fragments, which have withstood the 
disruptive effect of tidal torques when approaching the 
central star and have passed though the sink cell almost 
intact. These events are usually characterized by rather 
short rise and fall times due to the compact nature of the 
fragments triggering the bursts. 

The other type are closely-packed bursts, which occur 
one after another as shown in panels a, b, and f. As 
a rule, there is one primary, most energetic burst and 
a few secondary bursts of lesser amplitude. These clus¬ 
tered bursts are caused by infalling fragments that have 
started to disintegrate on their approach to the star due 
to the disruptive influence of tidal torques. An example 
of such a phenomenon is illustrated in Figure |T3| showing 
the gas surface density in model 1 in the inner region 
during a short time period corresponding to Figure |T2h- 
An extended fragment approaching the star is evident in 
the upper panels of the figure. At a distance of a few 
tens AU from the star, the fragment loses its roundish 
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Fig. 13. — Zoom-in onto a fragment approaching the central star. 
Shown is the gas surface density (in log g cm~^) duri ng a short 
time period corresponding to a clustered burst in Figure [T^ . The 
star is marked by the red circle in the coordinate center. 

shape, stretches into a clumpy filament, and hnally ac¬ 
cretes through the sink sell onto the star. It is interest¬ 
ing to note that another fragment seen in the bottom-left 
panel is hurled to a larger distance due to the gravita¬ 
tional exchange of angular momentum with the infalling 
fragment, a phenomenon that can potentially contribute 
to the survival probability of fragments in the disk. De¬ 
pending on the final structure of the filament, two or 
more bursts of varying amplitude can be triggered. 

The typical times between individual bursts in the 
clustered-burst event are 30-100 yr, too long to be firmly 
detected on the basis of available FUori observations. 
Nevertheless, the knotty structure seen in many proto- 
stellar jets and spacing between the knots are consistent 
with short periods of variability in the mass accretion 
rates as could be expected from clustered bursts. 

Another class of young pre-main-sequence stars, EX- 
ors (known after its prototype star EX Lupi), are known 
for exhibiting repetitive bu rsts of lesser magnitude on 
timescales of sever al years (|Hartmann &: Kenvonl [19^ 
lAiidard et al.ll201^ . Until recently, the repetitive nature 
of EXor bursts has been considered as a defining charac¬ 
teristic of such objects. The existence of clustered bursts, 
if confirmed, could further erode the ever shrinking gap 
between FUors and EXors. In this context, post-burst 
observations of known FUors using (sub)-millimeter in¬ 
terferometers may search for a clumpy filamentary struc¬ 
ture in the inner 20-30 AU caused by a disintegrated 
fragment approaching the star. A successful detection of 
such a structure would imply another outburst to come. 

7. THE FRACTION OF BURST-PRODUCING CORES 

Finally, based on our numerical simulations, we want 
to estimate the fraction of star-forming cores that is ex¬ 
pected to produce bursts after forming a star-disk sys¬ 
tem. Our model with core mass Mcore = 0.3 Mq does not 


undergo a burst mode driven by the formation of large 
fragments in the disk, however the models we run with 
Afcore = 1-1 AIq and greater mass do show significant 
bursts with mass accretion rate e xceed ing 10“"^ Mq yr“^. 
Figure 6 of iBasu fc VorobvovI (|2012t) shows a correla¬ 
tion between fragmentation events, core mass, and /?. 
Despite limitations of exploration of parameter space 
in these models, we can combi ne the information there 
with the observational finding of ICaselli et al.l (j2002il that 
P has a median value of 0.02 and most values in the 
range ~ 10“^ — 10“^ to conclude that most cores with 
Mcoie ^ 1.0 Mq will exhibit fragmentation and bursts 
and that most cores with Mcore < 0.5 Mq will not ex¬ 
hibit bursts. 

To make an estimate of what number fraction of star¬ 
forming events these represent, we can further employ 
the stellar initial mass function (IMF) of IChabrien (|2005ll 
and the idea that the core mass function is simply the 
same shape as the IMF but sca led up by a multip lica- 
tive factor of about 3 (see, e.g.. lAndre et al.l 1200911 . In 
this case, and using 0.075 M fn as a minimum stellar mass 
(|Chabrier fc Baraffell2000ll . we can identify the number 
of cores with Mcore > 1.0 Mq with the number of stars 
with mass m > 0.33 Mq. Their number fraction can be 
found eas ily using the cumu lative function of stars intro¬ 
duced bv iBasu et al.l (j2014ll in their equation (15): 




(13) 


-iexp ( ujno + 


■^erfc 


V2 


ln(m) - H o 
\/2ao 
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where erfc(x) = 1 — erf(a:), erf(a;) is the Error function, 
uj is the power-law index of the high mass tail of the 
mass distribution function (oc and fio and uo 

are additional parameters that help define the mean and 
variance of the distribution. 

Their best fit to the Chabrier IMF shows that about 
40% of all stellar objects have m > 0.33 Mq, hence 
presumed to originate from a core with mass Mcore > 
1.0 Mq. The criterion Mcore < 0.5 Mq for no bursts will 
translate to 0.075 Mq < m < 0.167 Mq in the stellar 
IMF, and this accounts for some 32% of all stellar ob¬ 
jects. If 40% of stellar objects are expected to definitely 
have bursts and 32% are not, this leaves an intermedi¬ 
ate 28% that may have bursts depending on the level of 
rotation in the initial core. Therefore, we can say that 
the fraction of star-forming cores that can be expected 
to display bursts after forming a star-disk system falls 
somewhere in the range of 40%-70%. 

8. CONCLUSIONS 

We have shown that variable protostellar accretion 
with episodic bursts is a robust property of protostel¬ 
lar collapse. This phenomenon is often present in a 
collapse environment in which a protostellar disk has 
a self-consistent mass loading from the core envelope. 
An accurate calculation of photospheric and accretion 
luminosities, as well as improved disk thermal physics 
and dust opacities, have refined our numerical hydrody¬ 
namics model, enabling better characterization of lumi¬ 
nosity bursts caused by disk gravitational fragmentation 
followed by fragments migrating onto the protostar. 
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In agreement with our previous studies, we found that 
an increase in the initial core mass and angular momen¬ 
tum favours the burst phenomenon, while higher levels 
of background radiation and magnetic fields can moder¬ 
ate the burst activity. A minimum mass infall rate onto 
the disk on the order of a few xlO“® Mq yr“^ is re¬ 
quired to generate the bursts. The main results can be 
summarized as follows. 

• A general correlation between the amplitude of 
time variations in M and the strength of non- 
axisymmetric perturbations in the disk, as defined 
by the global Fourier amplitudes, suggests a causal 
link between accretion variability and the disk 
gravitational instability. In our models, long-term 
variations in protostellar accretion are caused by 
the nonlinear interaction between different spiral 
modes in the gravitationally unstable disk, while 
episodic accretion bursts are triggered by fragments 
migrating onto the star. 

• Most luminosity bursts in our models occur in the 
partly embedded Class I phase, with a smaller frac¬ 
tion taking place in the deeply embedded Class 0 
phase and a few occasional bursts in the optically 
visible Class II phase. 

• The properties of the bursts are found to be in good 
agreement with those inferred for FU-Orionis-type 
objects (FUors). For instance, our models yield 
the mean luminosity and average burst duration of 
312 Lq and 39 yr, respectively, only a factor of 
1.5-2 higher than in FUors. 

• Depending on the ability of fragments to with¬ 
stand the tidal torques when approaching the star, 
two types of bursts can occur: the isolated and 
clustered ones. In the former, the fragment is 
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accreted almost intact, producing a well-defined 
burst, while in the latter the fragment is stretched 
into a knotty filament when approaching the star, 
thus producing a series of closely-packed bursts of 
varying amplitude separated from each other by 
several decades. 

• Adopting the stellar IMF of IChabri^ (|2005[ 1 and 
assuming that the core mass function is the same 
shape as the IMF but is scaled up by a factor 
of about 3, we estimate that about 40%-70% of 
the star-forming cores can be expected to exhibit 
bursts after forming a star-disk system. 

In the present work, the mass accretion rates M were 
calculated at the position of the inner sink cell at 6 AU. 
We have shown that the time behavior of M is similar 
at 8 AU and 12 AU, thus demonstrating that the accre¬ 
tion variability and the bursts are not the artifacts of the 
adopted inner boundary condition in our numerical sim¬ 
ulations. Nevertheless, more work is needed to develop 
a self-consistent link between the inner and outer disk 
regions. 

9. ACKNOWLEDGMENTS 

The authors are thankful to the anonymous referee 
for providing critical comments that helped to improve 
the manuscript. The authors are thankful to Isabelle 
Baraffe and Gilles Chabrier for providing the stellar evo¬ 
lution code. This work is partly supported by the RFBR 
grant 14-02-00719. The simulations were performed on 
the Shared Hierarchical Academic Research Computing 
Network (SHARCNET), on the Atlantic Computational 
Excellence Network (ACEnet), and on the Vienna Scien¬ 
tific Cluster (VSC-2). SB was supported by a Discovery 
Crant from the Natural Sciences and Engineering Re¬ 
search Council (NSERC) of Canada. 


Chen, H., Myers, P. C., Ladd, E. F., & Wood, D. O. S. 1995, ApJ, 
445, 377 

Chabrier, G. 2005, in The Initial Mass Function 50 years later, ed. 
E. Corbelli, F. Palla, and H. Zinnecker, (Springer, Dordrecht), 

41 

Chabrier, G. &; Baraffe, 1. 2000, ARA&A, 38, 337 
D’Antona, F., &; Mazzitelli, I. 1994, ApJS, 90, 467 
Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092 
Dunham M. M., Evans N. J., II, Terebey S., Dullemond C. P., 
Young C. H., 2010, ApJ, 710, 470 
Dunham, M. M., Wrobyov, E. 1. 2012, ApJ, 747, 52 
Dunham, M. M., et al. 2014, in Protostars and Planets VI, eds. 
H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 
University of Arizona Press, 195 

Evans, N. J., H, Dunham, M. M., Jorgensen, J. K., et al. 2009, 
ApJS, 181, 321 

Foster, P. N., Chevalier, R. A. 1993, ApJ, 416, 303 
Hartmann, L., & Kenyon, S. J. 1996, ARAA, 34, 207 
Johnson, B. M., Gammie C. F. 2003, ApJ, 597, 131 
Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 
1990, ApJ, 99, 869 

Kim, H. J., Evans, N. J., H, Dunham, M. M., Lee, J.-E., & 
Pontoppidan, K. M. 2012, ApJ, 758, 38 
Kospal, A., Abraham, P., Acosta-Pulido, J. A., et al. 2011, A&A, 
527, 133 

Larson, R. B. 1969, MNRAS, 149, 271 

Lee, J.-E. 2007, J. Korean Astron. Soc., 40, 83 

Machida, M. N., Inutsuka, S., Matsumoto, T. 2011, ApJ, 729, 

42 

Masunaga, H., & Inutsuka, S.-I. 2000, ApJ, 531, 350 






18 


Vorobyov & Basu 


Meru, F., & Bate, M. R. 2012, MNRAS, 427, 2012 
Nakano, T., Nakamura, T. 1978, PASJ, 30, 671 
Nayakshin, S., & Lodato, G. 2012, MNRAS, 426, 70 
Ohtani, T., Kimura, S. S., Tsuribe, T., Vorobyov, E. I. 2014, PASJ, 
in press 

Padoan, P., Haugbolle, T., Nordlung A. 2014, arXiv: 1407.1452 
Pension, M. V. 1969, MNRAS, 144, 425 

Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., Si 
Denzmore, P. 2006, ApJS, 167, 256 
Quanz S. P., et al. 2007, ApJ, 668, 359. 

Rice, W. K. M., Mayo, J. H., Si Armitage, P. J. 2010, MNRAS, 
402, 1740 

Scholz, A., Froebrich, D., Si Wood, K. 2013, MNRAS, 430, 2910 
Michael, S. Si Durisen, R. H. 2010, MNRAS, 406, 273 
Millan-Gabet R., et al. 2006, ApJ, 641, 547 

Semenov, D., Henning, Th., Helling, Ch., Ilgner, M., Si Sedlmayr, 
E. 2003, ASiA, 410, 611. 

Scholz A., Froebrich, D., Wood, K. 2013, MNRAS, 430, 2910 
Shu, F. S. 1977, ApJ, 214, 488 

Stamatellos, D., Whitworth, A. P., Si Hubber, D. A. 2011, ApJ, 
730, 32 


Tsukamoto, Y., Takahashi, S. Z., Machida, M. N., Inutsuka, S. 
2014, ArXiv:1404:7271 

Visser, R., Si Bergin, E. A. 2012, ApJ, 754, 18 
Vorobyov, E. 1. 2009, ApJ, 704, 715 
Vorobyov, E. L 2010, ApJ, 723, 1294 
Vorobyov, E. L 2011, ApJ, 729, 146 
Vorobyov, E. L 2013, ASiA, 552, 129 
Vorobyov, E. L, Si Basu, S. 2005, MNRAS, 360, 675 
Vorobyov, E. L, Basu, S., 2006, ApJ, 650, 956 
Vorobyov, E. L, Si Basu, S. 2009, MNRAS, 393, 822 
Vorobyov, E. L, Si Basu, S. 2010a, ApJ, 719, 1896 
Vorobyov, E. L, Si Basu, S. 2010b, ApJL, 714, 133 
Vorobyov, E. L, Baraffe, L, Harries, T., Chabrier, G. 2013, A&A, 
557, 35 

Vorobyov, E. I., Zakhozhay, O. V., Dunham, M. M. 2013, MNRAS, 
433, 3256 

Zhu, Z., Hartmann, L., Si Gammie, C. 2009, ApJ, 694, 1045 
Zhu, Z., Hartmann, L., Nelson, R. R, Si Gammie, C. F. 2012, ApJ, 
746, 110 


